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ABSTRACT 

We consider the inner regions of accretion disks surrounding black holes and neutron 
stars and investigate the nonlinear time dependent evolution of thermal- viscous instabilities. 
The viscous stress is assumed to be proportional to the gas pressure with the viscosity 
parameter formulated as a = mm[ao{h/r) n , a max ], where h is the local scale height, r is the 
distance from the central compact object, and n, «o and a max are constants. It is found that 
the disk is unstable for a sufficiently sensitive to h (n ^ 1.2). The instabilities are globally 
coherent in the entire unstable region of the disk and, depending on the viscosity parameters, 
the time variability of the mass accretion rates are manifested as periodic or quasi-periodic 
oscillations. We show that, the low frequency (~ 0.04 Hz) quasi-periodic oscillations (QPOs) 
discovered recently in some of the black hole candidates (Cyg X-l and GRO J0422+32) and 
a low mass X-ray binary (Rapid Burster MXB 1730-335) may be explicable by the thermal- 
viscous instabilities in accretion disks. The observations of QPOs place constraints on the 
viscosity parameters and suggest that (n, a ) ~ (1.6, 30) for the Rapid Burster with a 1.4 M 
neutron star. In the case of black hole candidates, the dependence of a on h/r is less steep 
corresponding to n ~ 1.2 — 1.3 for black holes less than 10 M . 

Subject headings: 
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1. INTRODUCTION 

The quasi-periodic oscillation (QPO) phenomenon in low mass X-ray binaries (LMXBs) 
and galactic black hole candidates (GBHCs) provides a powerful tool to probe the temporal 
behavior of these systems. An understanding of the observed luminosity variations may 
place important constraints on the physics of accretion in these objects. QPO phenomenon 
in bright LMXBs have been known since the mid 1980s (see Lewin, van Paradijs, & van 
der Klis 1988; van der Klis 1989). Two distinct kinds of QPOs have been classified, one is 
the intensity-dependent QPOs characterized by frequencies ~ 20 — 50 Hz, and the other is 
the intensity-independent QPOs with frequencies of ~ 5 — 10 Hz. The former QPOs are 
called horizontal-branch QPOs whereas the latter QPOs are called normal branch QPOs. 
The horizontal-branch QPO phenomenon has been interpreted in terms of an oscillatory 
variation in the luminosity due to modulations in the mass accretion rate and the QPO 
frequencies are identified with a beat frequency corresponding to the difference between the 
Keplerian orbital frequency at the magnetopause and the rotation frequency of the central 
neutron star (Alpar & Shaham 1985; Lamb, Shibazaki, Alpar, & Shaham 1985). The normal 
branch QPOs, on the other hand, have been modeled as a radiation instability in a spherical 
accretion flow (Fortner, Lamb, & Miller 1989) or as a variation in the vertical structure of an 
accretion disk (Alpar, Hasinger, Shaham, & Yancopoulos 1992) for accretion rates near the 
Eddington limit. In contrast to the model for the horizontal branch QPOs, the normal branch 
QPOs are interpreted in terms of optical depth variations of a nearly constant luminosity. 

More recently, another kind of QPOs from X-ray binary sources have been discovered. 
They are characterized by low frequencies of ~ 0.04 Hz in both neutron star and black hole 
candidate systems. For example, such QPOs have been found in the hump phase of the 
persistent emission of the Rapid Burster MXB 1730-335 (Lubin et al. 1992, 1993), in the 
low state of Cyg X-l (Angelini & White 1992; Vikhlinin et al. 1992a, 1994; Kouveliotou et 
al. 1992a), and in GRO J0422+32 (Vikhlinin et al. 1992b, Kouveliotou et al. 1992b, Pietsch 
et al. 1993). While the Rapid Burster is believed to be a neutron star, the other two sources 
are black hole candidates. If a black hole is involved, the mechanism deemed responsible for 
the production of these low frequency QPOs may be different from those proposed for the 
other classes of QPOs (the above models for the horizontal or normal branch QPOs require 
the central objects to be neutron stars) since such oscillations would then be independent 
of the nature of the compact object. Thus, it is natural to explore the possibility that 
disk instabilities are involved in the phenomenon. Among the various models proposed 
for the general QPO phenomenon is a model suggested by Abramowicz, Szuszkiewicz, & 
Wallinder (1989) who proposed that thermal and viscous instabilities in accretion disks 
could be important for QPOs with frequencies of ~ 1 Hz for LMXBs and ~ 0.1 Hz for 
GBHCs. 

It is well known that the inner regions of viscous accretion disks surrounding black holes 
or neutron stars may suffer thermal- viscous instabilities when radiation pressure is important 
(see e.g. , Piran 1978). For a standard a-model of accretion disks (see Shakura & Sunyaev 
1973 ) with the viscous stress, r, proportional to the total pressure, p, i.e. , r = — ap, 
where a is a constant, disks are thermally and viscously unstable if the pressue is radiation 
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dominated (Lightman & Eardley 1974; Shakura & Sunyaev 1976). However, if the viscous 
stress is proportional to the gas pressure, p g , only, i.e. , r = —ap g , and a is a constant (see, 
e.g. , Lightman & Eardley 1974; Coroniti 1981; Stella & Rosner 1984), then the disk is always 
secularly stable. 

Time dependent calculations of the thermal-viscous instabilities have been carried out to 
examine the global behavior of the disk (see, e.g. , Taam & Lin 1984; Honma, Matsumoto & 
Kato 1991; Lasota & Pelat 1991). It was revealed that, for the r = — ap prescription, the 
disk is globally unstable and large-amplitude, burst-like fluctuations of the disk luminosity 
arise if the mass accretion rates are high enough that the pressure in the disk is radiation 
dominated. However, this strong instability may be reduced or even suppressed by modifying 
the a-viscosity prescription. For example, with r = —ap/3 m (see Abramowicz, Szuszkiewicz 
& Wallinder 1989; Szuszkiewicz 1990), where j3 = p g /p, the disk is stable for m = 0.5 (Taam 
& Lin 1984; Honma et al. 1991). For m = 0.25, the amplitude and the time scale of the 
fluctuations of the disk luminosity is largely reduced (Honma et al. 1991). The stabilization 
in the case of m — 0.5 or effectively r = —a^/pp g in Keplerian disks, is not predicted by 
local linear analysis. It is due to the advection of energy by material motion in the radial 
direction (Taam & Lin 1984). 

Notwithstanding the lack of a fundamental theory for the viscosity in accretion disks, 
the two phenomenological viscosity prescriptions, i.e. , r = — ap and r = — ap g , have some 
theoretical support based on a turbulent viscosity picture (see Shakura & Sunyaev 1973) and 
magnetic field induced viscosity picture (see, e.g. , Lightman & Eardley 1974; Coroniti 1981; 
Stella & Rosner 1984), respectively. On the other hand, the disk instabilities based upon a 
viscous stress proportional to the total pressure may be too strong to be relevant to the low 
frequency QPOs, and disks with the gas pressure prescription are stable and can not intro- 
duce any variability within the assumed framework. A possible resolution to this dilemma is 
to introduce a mixed prescription. In other words, as mentioned above, one may introduce a 
dimensionless parameter m between and 1. However, there is another possibility, namely, 
that the viscous parameter, a, is not a constant. In fact, a non-constant a is usually required 
in the studies of accretion disks in dwarf novae. In that case, the disk temperatures are low 
and radiation pressue is unimportant in comparision with the gas pressure. However, due to 
the sensitivity of the opacity on the temperature, a thermal instability occurs. These insta- 
bilities have been applied to explain the dwarf nova phenomena, and in order to understand 
both the recurrence time and the duration of the outbursts, different values of the viscosity 
parameter a are required during the two states. One widely applied formula which can 
approximate such a variation is the form of a — a (h/r) n (e.g. , Meyer & Meyer-Hofmeister 
1983; Duschl & Livio 1989), where h is the local scale height of the disk and r is the distance 
from the compact object. This form also has a theoretical basis from scaling arguments for 
magnetic field induced viscosity. In particular, for a mean helicity determined by cyclonic 
convective motions, Meyer & Meyer-Hofmeister (1983) find n = 1.5, whereas for a helicity 
due to internal waves, Vishniac & Diamond (1992) conclude that n — 4/3. In addition, it 
was shown that, the effective a due to spiral shock waves varies like (h/r) 1 - 5 , even though 
its magnitude is very small (see Spruit 1987). 
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The {h/r) n dependence has a destabilizing effect (assuming n is positive) on the disk if 
it is applied to the high temperature inner regions. Based on the previous nonlinear studies, 
the r = — ap model is excluded. We apply this form for a to the most stable model, i.e. , 
r = —ap g . In this case, it has been shown that the disk is locally unstable (Meyer 1986) for 
n ^ 0.75 in the limit that radiation pressure in dominant. In a recent study, Milsom, Chen, 
& Taam (1994, hereafter MCT) generalized the results of Meyer (1986) and, furthermore, 
pointed out that, in this model, unlike the r = —ap prescription, the instability may take a 
milder form and non burst-like oscillations may exist in some circumstances, which may be 
applied to an explanation of the low frequency QPOs. 

In this paper, we investigate the global evolution of the thermal-viscous instabilities in 
accretion disks in a time dependent approach to determine their viability as a mechanism 
for low frequency QPOs. It is the purpose of this paper to report on the detailed nonlin- 
ear time dependent calculations of accretion disks which suffer thermal-viscous instabilities 
based upon the above modified form of the viscosity prescription. In the next section, the 
fundamental equations and the approach we adopt are outlined. In §3 a brief review of the 
local thermal and viscous instability in accretion disks is presented. The detailed numerical 
results of the global evolution are given in §4, and the implications of our results and their 
possible applications for interpretation of observations will be discussed in the last section. 

2. FORMULATION 

We focus only on the thermal-viscous instability which may arise in accretion disks. The 
axisymmetric inertial-acoustic instabilities, which have a much shorter time scale and may 
exist possibly in the innermost regions of the disk (see, e.g. , Matsumoto, Kato, & Honma 
1989; Nowak & Wagoner 1991; Wallinder 1991; Chen & Taam 1993), are not considered 
here. Non-axisymmetric or self-gravitating effects (see, e.g. , Papaloizou & Pringle 1984, 
1985, 1987), are also neglected. In other words, we assume a Keplerian disk, which is 
axisymmetric, non self-gravitating, optically thick and geometrical thin so that it can be 
described by the vertically integrated equations. In this approximation, the surface density 
of the disk, E, at a given cylindrical radius, r, can be obtained by combining the conservation 
equations of mass and angular momentum. It is given by the standard time dependent mass 
diffusion equation (Lynden-Bell & Pringle 1974) as 

d^_3d_ 
dt r dr 

where t is the time, and v is the effective kinematic viscosity. In equation (1), the Newtonian 

Keplerian angular velocity, Q = ^jGM/r 3 , is implicitly assumed, where M is the mass of the 
central compact object. Relativistic effects on the gravitational potential are important near 
the inner edge of the disk, but we anticipate that the long time scale (~ 20 s) associated 
with the low frequency QPOs dictates that the region of instability lies at radii sufficiently 
large compared to the innermost disk that the results we obtain will be insensitive to the 
neglect of these effects. 
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The mid-plane temperature, T, of the disk is governed by energy conservation which is 
expressed as 



cy:r\( dlnT + v dlnT \ (r i)P lnS + v ainE 91nir 
C ^ T [[-W + Vr ^~ J" (r3 " 1} [~W + Vr —r dT / 

where V r , T 3 , and C„ are the radial velocity, the adiabatic exponent and the specific heat at 
constant volume respectively. H is the half-thickness of the disk and is defined as 

The terms on the right hand side of equation (2) describe the local heating and cooling 
processes, where F + is the viscous dissipation rate, F~ is the cooling rate in the vertical 
direction, and F r is the radiative energy flux in the radial direction. The heating rate and 
viscosity is related by 

And the energy transport flux in the radial direction is written as 

F r = -*HF- d -^. (5) 

The effective kinematic viscosity is parameterized in terms of the a model and is assumed 
to be 

2 

v = -ac s h(p g /p), (6) 

where c s and h are the local sound speed and scale height respectively. The viscosity param- 
eter, a, is formulated as 

a = mm[a (h/r) n , a max ], (7) 

where a is a constant and a max (< 1) represents the effect of saturation of the anomalous 
turbulence. 

In the radiative diffusion and vertically average approximation, the viscosity, P, is cal- 
culated by using the half-thickness of the disk, H, as the scale height h and y / 2~^p/S as 
the sound speed c s . This form of viscosity leads to a viscous stress corresponding to the 
r = —ap g model in the Keplerian approximation, and it has a corresponding heating rate of 
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In this approximation, the cooling rate can be expressed as 



F~ = 



4acT 4 
3kS 



(9) 



where k is the sum of electron scattering and free-free opacities. 

The heating and cooling rates in the energy equations can also be calculated directly 
from the detailed vertical structures of the disk. In practice in the time dependent study of 
dwarf nova accretion disk models, the above simplified heating and cooling rates in equations 
(8) and (9) are usually not applied directly because they differ from the vertical structure 
results qualitatively due to the temperature sensitivity of the opacity. Many schemes have 
been introduced to deal with this situation (see, e.g. , Papaloizou, Faulkner & Lin 1983; Mi- 
neshige & Osaki 1983; Meyer & Meyer-Hofmeister 1984; Cannizzo, Wheeler & Polidan 1986; 
Mineshige 1986). For example, Faulkner, Lin & Papaloizou (1983) and Papaloizou, Faulkner 
& Lin (1983) used a modified approximate expression for the cooling rate. Mineshige & 
Osaki (1983) tabulated both the cooling and heating rates from detailed vertical structures. 
To construct the disk vertical structures, four equations (which are, in the steady state ap- 
proximation, hydrostatic equilibrium, energy transport and energy conservation equations) 
are required to obtain the solutions of density, temperature and vertical energy flux with re- 
spect to the height from the disk mid-plane. In the approach of Mineshige and Osaki (1983), 
the condition of thermal equilibrium in the vertical direction (i.e. , the energy conservation 
equation) is relaxed, and instead, a relationship between the energy flux and the height is 
assumed. 

For disks surrounding black holes and neutron stars, MCT constructed detailed steady 
state vertical structures including convection and showed that convection has a stabilizing 
effect on the disk and changes the disk structure significantly. Therefore, we choose to use the 
cooling and heating rates from the detailed vertical structure models instead of the simplified 
rates above. To tabulate the cooling and heating rates, a different approach from that of 
Mineshige & Osaki (1983) is adopted as follows. 

We assume that the departure from thermal equilibrium in the vertical direction has a 
form of 



where, 5 is a dimensionless constant and is introduced to represent the effect of advection of 
energy by material motion and radiative energy transport in the radial direction. For 8 > 1, 
the radial energy transport is a heating process, and for 5 < 1, it is a cooling process. The 
energy equation in the vertical direction becomes 



where z is the height from the disk mid-plane, F, p and e„ are the total flux, the density and 
the viscous energy generation rate at height z respectively. Note that, in the calculation of 
e u , which is 9/4z/f2 2 , the viscosity is computed with the local adiabatic sound speed and a 
modified local scale height and they are functions of z (see MCT). Now, if we identify the 



F~ = 5F+, 



(10) 




(11) 
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temperature T(z) at the mid-plane as T c , the total flux F(z) at the surface as F s and the 
column density E(z) at the surface as E s , then we have the following relations for E, T, F~ 
and F+: 

S = 2E„ 

(12) 



T = T C , 



F- = 2F S , 
F + = 2F s /5. 

In general, 5 is unknown and F~ and F + can not be obtained from E and T. However, 
for a given viscosity prescription and specified radius, it is found that two relations for F~ 
and F + exist which are approximately independent of 5. These relations can be formally 
expressed as 

f / 1 (lnF-) = lnF- 

\ / 2 (lnF-) = F+/F + - 1 ' 

As an example, f\ and f 2 are shown in Figure 1(a) and (b) respectively for a specified set 
of model parameters and radius. In each plot, there are three curves of 8 — 0.8, 1.0 and 
1.5 (corresponding to dotted, solid and dashed lines respectively), however, the difference 
between them is very small. The f 2 curve is Z-shaped, and is related to the S-shaped 
relation between the temperature and surface density of the disk which will be discussed in 
the following sections. 

These two functions are tabulated with respect to F~ at specified radii for a given vis- 
cosity prescription. For the time dependent evolution, F~ and F + are calculated from the 
two variables E and T (see eqs. [8] and [9]), and linear interpolation is used to obtain fi 
and f 2 from the tables to determine the modified values of F~ and F + , i.e. , 

(14) 

The time dependent equations (1) and (2) are solved to calculate E and T via an explicit 
method. The accretion disk is divided into 41 grid points distributed equally on logarithmic 
scale ranging from an inner boundary at Ar g to an outer boundary at 300 r g , where r g is 
the Schwarzschild radius (r g = 2GM/c 2 ). The initial structure of the disk is given by the 
thermal equilibrium solutions (5 = 1) of the detailed vertical structures (see MCT). At every 
grid radius a series of models of steady state vertical structures with different mass accretion 
rates are constructed, and at the same time, functions of f\ and f 2 are tabulated. For the 
initial values of E and T with the same steady state mass accretion rate at each grid point, 
linear interpolation from the vertical structure models of the same grid radius is applied. 
The boundary conditions are chosen to correspond to zero gradients at the innermost radius 
and to a fixed temperature and surface density at the outermost radius. 

The disk model is determined by the mass of the central object, M, the mass accretion 
rate, M, and the viscosity parameters, a , n, and ct max . We adopt the mass as 1.4 M Q for 
neutron stars and 1OM for black holes. The mass accretion rate is measured in units of 
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the Eddington limit defined as Mgdd — 4:7iGM/(K, e ce), where n e is the electron scattering 
opacity and e is the efficiency for the conversion of rest mass energy into radiation taken to 
be 1/6 for a neutron star and 1/16 for a black hole. 

3. THERMAL AND VISCOUS INSTABILITIES 

Thermal instability occurs when the heating process is more efficient than the cooling 
process. In an optically thick and geometrical thin viscous accretion disk, if we consider only 
the local processes, then the heating is due to the viscous dissipation and the cooling is due 
to the radiative diffusion (the convection is neglected). The general criterion for thermal 
instability may be expressed as (see Piran 1978) 



<91nF+ 



d\nH 



dlnF~ 

> 



d\nH 



(15) 



The local thermal instability of accretion disks for which the viscous stress is given in 
the form of r = —atQ{h/r) n p g has been examined by Meyer (1986) in the limit in which 
radiation pressure dominates gas pressure and it was shown that thermal instability occurs 
for n > 0.75. In a recent study, MCT generalized the analysis to arbitrary values of the 
ratio of gas pressure to total pressure (($) and derived the following condition for thermal 
instability: 

"<^)=^- < 16 > 

Furthermore, it was shown that this condition is identical to the viscous instability condition 
which can be inferred from the slope of the relation between the mass accretion rate and 
surface density at a fixed radius. That is, for ^ > the disk is locally stable whereas for 

< it is locally unstable (see Bath & Pringle 1982). In particular, in the polytropic 
approximation, the slope is given as (see MCT) 

dlnM 5 + (5 + n)(3 

dlnS ~ 4n-3-3(n+l)/3' ^ ' 

From equation (17), for a given n (less than 6), it is seen that, in the parameter plane of mass 
accretion rate and surface density, the slope of M(E) relation is positive if the gas pressure 
is much larger than the radiation pressure, and which corresponds to the stable regime. At 
the point (3 = (3 C , the slope changes sign from ^ > to ^ < 0, and the condition of 
j3 < p c corresponds to the unstable regime. 

In the detailed vertical structures of accretion disks, the above feature is qualitatively 
recovered. However, another turning point was revealed, which is due to the saturation 
of the a parameter (see MCT), and, hence, an S-shaped curve is formed in the plane of 
mass accretion rate and surface density. The upper branch is stable because in that regime 
(llnir — 5/3) the disk is approximately described by the r = — a max p g prescription, and no 
dependence of the viscous parameter on the local scale height exists. 
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An S-shaped relation of temperature with respect to surface density follows from the 
relation between T and M, specifically, 

M oc F~ oc T 4 /S. (18) 

It is also seen (eq. [17]) that, for larger n, \^-\ or |^| at the middle branch of the S-curve 
is smaller and so the S-curve becomes steeper, which results in a greater tendency toward 
instability. This can also be explained through the different sensitivities of the cooling and 
heating rates on the temperature. For the the viscous heating rate, we have, approximately, 



<91nF+ 



<91nT 

which becomes larger for larger n. On the other hand, 



,4-3/3 7-7/5 

^ (2 (19) 



d\nF~ 



<91nT 



4, (20) 



so the temperature sensitivity of the cooling rate is approximately a constant. 

A local thermal- viscous instability at a fixed radius follows from the evolutionary trajec- 
tory in the (T, E) plane. For example, as the mass accretion rate increases while the disk is on 
the lower stable branch, the surface density increases until the turning point where (5 = f3 c is 
reached. The subsequent evolution is expected to lead to heating of the disk until the upper 
stable branch is reached. The mass flow rate corresponding to this state, however, is higher 
than that corresponding to the lower stable branch and therefore S decreases. Eventually, 
the evolutionary path reaches the upper turning point and the disk undergoes a transition 
to the lower stable branch, where the mass flow rate is less than the assumed steady state 
rate. As a result, matter accumulates, thereby increasing X and the cycle begins anew. To 
determine whether the disk follows the steady state disk curves, as outlined above, we turn 
to the description of the global evolution. 

4. GLOBAL EVOLUTION 

The parameters characterizing each of the model sequences as well as the results of the 
simulations are summarized in Table 1. Here, M is in units of the Eddington value. For 
the unstable models, the frequency of the oscillation, /, and the relative amplitude of the 
luminosity fluctuation, AL/L, are also listed. 

4.1 Neutron Stars 

The existence of an S-shaped relation between T and £ at a specified radius indicates 
thermal and viscous instabilities whenever the steady state solutions of T and E are located at 
the middle branch. However, these instabilities are local and may be stabilized by nonlinear 
effects especially if the instability is restricted to only a small region of the disk. Here, we 
present the results of the time dependent calculation of model sequence 1 (n = 1.1) with a 
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mass accretion rate of 0.3 Mgdci and a central object of mass equal to 1.4 M to demonstrate 
the global stabilizing effects. The steady state (i.e. , the initial state) is indicated by open 
circles in the T and E plane at four specified radii (see Fig. 2). The ratios of gas to total 
pressure, (3, are also listed for each radius. In the radiative diffusion approximation, for 
n = 1.1, the critical value of the ratio is, f3 c = 0.2121 (see eq. [16]). Thus, the disk is 
predicted to be unstable in three of the four radii displayed. However, it is seen that, only 
two open circles in the inner radii are on the unstable middle branches, which confirms 
the results of MCT that the convection has stabilizing effects and a greater contribution of 
radiation pressure is required to make the disk locally unstable. The nonlinear calculation 
shows further that the global effects, i.e. , the advection of energy and radiation transport 
in radial direction, stabilize the entire disk, even though it is locally unstable in the inner 
regions. As is shown in Figure 2, the evolution models are seen to be approximately the 
same as that of the initial steady state (the last evolution model is indicated by the cross). 
The variation of disk luminosity, calculated by integrating the flux F~ over the area of the 
disk surface, with respect to time is displayed on Figure 3. Obviously, it evolves to a steady 
state after a short initial transient stage of less than 10 seconds. Additional calculations 
have been performed for different mass accretion rates ranging from 0.01 to 1 times the 
Eddington value, and the disk is always found to evolve to a stable steady state structure. 
These results do not confirm local analysis which predicts instability for n ^ 0.75 and reveal 
that stabilization for n ^ 1.1 results from the global effects of non-local energy transport 
(see also Taam & Lin 1984). 

The behavior of the accretion disk is a sensitive function of the parameters characterizing 
the viscosity prescription (see §3). To strengthen the disk instability, we increase the viscosity 
parameter index n to 1.3 and keep the other viscosity parameters fixed as that of model 
sequence 1. We denote that model as sequence 2. The results for M = 0.27 MEdd are shown 
in Figure 4 and Figure 5a. It is seen that, the evolution paths of T and £ at radii 5.53 r g 
and 11.77 r g ( see Figs. 4a, b) are both very narrow loops moving in a counterclockwise 
direction. The narrow loops represent a weak instability, which can also be seen from the 
small amplitude periodic oscillations of the disk luminosity (Fig. 5a). In particular, the 
frequency of the oscillation is about 0.143 Hz and the amplitude is about 22%. At the other 
two larger radii, 31.10 r g and 59.43 r g , the steady state solutions of T and E are located at 
the lower branch of the S-curve (see Fig. 4c, d), and they are stable. The instability becomes 
stronger at higher mass accretion rates (M = 0.3MEdd) because radiation pressure becomes 
more important and the unstable regions of the disk is wider. Therefore, a lower frequency, 
larger amplitude oscillation results (see Fig. 5b). On the other hand, for mass accretion rates 
less than 0.21 MEdd, the disk is eventually stablized (Fig. 5c). 

The low frequency QPOs of about 0.04 Hz from the Rapid Burster imply a relatively small 
effective a parameter. To obtain a smaller effective a parameter, one can either increase n or 
decrease cto- However, for a decrease in ao, the surface density increases and, for a given M, 
radiation pressure becomes less important in the disk. Instability can only, therefore, occur 
at higher mass accretion rates. Since the persistent mass accretion rate in the Rapid Burster 
is probably low ( ^ 0.2 MEdd), a larger n is suggested. Furthermore, the parameter a ma x 



10 



determines the range of mass accretion rates over which the disk is unstable. For smaller 
a max , the range of mass accretion rates is smaller. The absence of low frequency QPOs 
during the type II X-ray bursts may require a smaller a max . 

In model sequences 5 — 8, n — 1.6, ao — 30 and a max = 0.2 are applied. The larger 
value of n makes the disk more locally unstable and also results in a wider unstable region 
of the disk. This effect, in addition to a smaller effective viscosity, leads to a longer viscous 
time scale. We first examine sequence 5 with M = 0.14M Edd . Figure 6 shows the evolution 
path in the temperature and surface density space at four specified radii. As expected, at 
the inner radii, the path has larger loops because of a stronger instability. However, the 
evolution path does not jump vertically from the first turning point to the upper branch. 
Instead, the trajectory takes an intermediate route between the upper and middle branches 
(closer to the latter one). The evolution path returns to the lower branch by following a 
route also closer to the middle branch. On the other hand, it does follow the lower branch 
approximately because the effects of radial advection are negligible there. At a larger radius, 
the loop is smaller, and eventually, it dissapears. It should be noticed that, in this case, the 
unstable region is wider than indicated by the local instability analysis. As seen in Figure 6c, 
at radius 31.10 r g , the initial model is located at the lower stable branch, however, due to the 
global radial motions, perturbations in the inner regions move outwards beyond this point 
and the stable region becomes unstable. In this model sequence, the unstable region of the 
disk is spatially confined inside 40 r g . 

The light curve resulting from the mass flow modulations in the accretion disk is illus- 
trated in Figure 7a. The period of the oscillations is found to be ~ 21.5 s and the amplitude 
of the fluctuations correspond to ^ ~ 67%. 

To determine the sensitivity of the results to the mass accretion rate, the mass accretion 
rate was decreased to M/M Edd = 0.12 in sequence 6. The disk is unstable as evidenced by the 
luminosity variations illustrated in Figure 7b. The fluctuations are manifested as oscillations 
and are similar to Figure 7a. It is found that the period has decreased to lis, and the relative 
luminosity amplitude has decreased, to ~ 20%. However, the relation that the frequency of 
the oscillation decreases and the relative luminosity fluctuation amplitude increases as the 
mass accretion rate is increased, holds only over a limited range in mass accretion rates for 
which the disk is unstable. For higher mass accretion rates, the amplitude of the oscillations 
may decrease and eventually the disk becomes stable due to the saturation of the viscosity 
parameter, a. For a set of viscosity parameters (n, ao, a max ) ~ (1.6,30,0.2), it is found 
that the accretion disk is stable for M ^ 0.11M Edd (sequence 8) and for M ^ 0.47M Edd 
(sequence 9). 

For a large viscosity parameter index, such as n = 1.6, the variation of the disk luminosity 
may exhibit different evolution patterns in some range of mass accretion rates. One example 
is shown in Figure 7c where two local maxima can be seen. 

4.2 Black Hole Candidates 

To determine the dependence of the results on the mass of the compact object we in- 
creased M to 10 M & to model the evolution of an accretion disk surrounding a black hole. 
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In order to produce oscillations at frequencies ~ 0.04 Hz, the effective viscosity parameter, 
a, must be larger than that of the neutron star case since the absolute size of the inner disk 
is larger. We first keep ao = 30 unchanged, but use n = 1.2 (see model sequences 10-13, 
Table 1). The instability is moderately mild, as shown by the disk luminosity oscillations in 
Figure 8a. For example, in the case of sequence 10 with mass accretion rate of 0.076 M^dd, 
the oscillation has a frequency of 0.026 Hz, and a small amplitude of AL/L ~ 0.083. For 
an increase of mass accretion rate, the frequency of the oscillation becomes even lower (se- 
quence 11). The range of mass accretion rates in which the disk is unstable is about a factor 
of 3, specifically, 0.07 ^ M/M Edd ^ 0.20 (sequences 12 and 13). In order to obtain a higher 
frequency, i.e. , 0.04 Hz, ao = 50 is applied in sequence 14-16. The frequency of the luminos- 
ity oscillations of sequence 14 with mass accretion rates 0.076 M Edd is increased to 0.035 Hz 
(see Fig. 8b). However, the range of mass accretion rates for which instability is indicated is 
very narrow in comparision with the model sequences with a = 30 (i.e. , sequences 12-13), 
specifically, 0.071 ^ M/ME dd ^ 0.084. The smaller upper limit on the mass accretion rate 
is due to the larger ao which makes the viscosity parameter saturate more easily. 

The mild instability indicated by n = 1.2 suggested that this value of the viscosity 
parameter index is near the lower limit for an unstable disk. This is confirmed by model 
sequence 17, which is performed to shown that n = 1.1 is stable for a 1OM black hole in a 
wide range of mass accretion rates, specifically, 0.005 ^ M/M Edd ^ 2.0. For a larger n, the 
strength of the instability become stronger and for a fixed ao the period of the oscillation is 
longer. For a 10 M black hole, the possible modulation frequency of oscillations of the disk 
luminosity may always be lower than the 0.04 Hz QPOs observed. Thus, we repeated the 
simulations of the models with viscosity parameters of n = 1.2, ao = 30, and a max = 1.0, 
but with a 5 M black hole. In this case, the range of mass accretion rates in which the disk 
is unstable is narrower than in that of a 10 M black hole (see sequences 12-13), specifically, 
0.085 ^ M/M Edd ^ 0.18. The frequency of the oscillations, as expected, is higher. An 
example is given in model sequence 18 with M = 0.1M Edd , where, the frequency is about 
0.05 Hz. 

In the case of a 5 M black hole, a wider range of mass accretion rates in which the 
disk is unstable can be obtained with n = 1.3 (sequence 19-21), and the time scale of the 
oscillation is in the required range, i.e. , ~ 0.04 Hz (see Fig. 8c for sequence 19) 

5. DISCUSSION 

We have demonstrated by global analysis that accretion in a geometrically thin, optically 
thick disk surrounding either a neutron star or a black hole can be unstable. The instability is 
thermal in origin and results from the inability of the disk to maintain a local thermal balance 
whenever radiation significantly contributes to the pressure in the disk. The strength of the 
instability is determined by the sensitivity of the viscous heating rate to the temperature, and 
is greater for larger n. The instability may be restricted to a relatively narrow spatial extent 
in the disk, and the nonlinear calculations reveal that the instability results in luminosity 
oscillations rather than bursts. Such oscillations may exhibit a quasi-periodic behavior if the 
mass flow rate entering the inner region of the disk is not strictly constant. 



12 



The thermal-viscous instability of the accretion disks can be understood through the 
steady state T ~ S relation curve (or effectively, the relation between M and S). The S- 
shaped relation for T(S) exists in slim accretion disk models (see e.g. , Abramowicz, Czerny, 
Lasota & Szuszkiewicz 1988; Szuszkiewicz 1990; Chen & Taam 1993), with its upper stable 
branch due to radial advection, which is a cooling process. The stabilization occurs only at 
super-Eddington rates. The S-shaped relation of T and E is also well known in disk models 
of dwarf novae. In that case, the stable upper branch represents a nearly fully ionized 
radiative disk structure with the gas pressure dominant. The middle unstable branch is due 
to the partial ionization of the material, which results in a sensitive temperature dependent 
opacity. The lower stable branch corresponds to cool non-ionized disk structure. It has 
been commonly assumed that the evolutionary path in T and S follows the stable lower and 
upper branches but does not follow the unstable middle branch. Instead, the path makes 
a vertical transition at constant £ from the lower turning point to the upper branch, and 
a downward transition from the upper turning point to the lower branch. This assumption 
is usually confirmed by the global time dependent calculations (see review by Cannizzo, 
1993). However, it is not obvious that is always the case because the evolution path may not 
neccessarily follow the steady state curve, which is a result of the vertical thermal equilibrium. 
For example, the evolution path will be affected by the radial motion of the disk which could 
even stablize it. This has been demonstrated by our results and, in fact, it has been also 
shown by the calculations of Honma et al. (1991) in the slim disk approximation. Our results 
show that, disks with viscosity parameter index n & 1.1 are globally stable. Our results also 
reveal that the evolutionary paths in the T and E plane are loops and the size of the loop 
reflects the strength of the global instability of the disk. 

Even though the time scales of the thermal and viscous instability at different radii of 
the disk are different, the global evolutions of the disk show that the evolution time scale is 
not determined locally and the instabilities are globally coherent. To see this, we plot the 
time variations of the radial distribution of the mass accretion rate of model sequence 19 (see 
Table 1) in Figure 9. The instability is initiated in the inner region where the contribution of 
radiation to the total pressure is the greatest and, as a result, a large amount of mass moves 
inwards at a high radial velocity. This results in a low surface density and a transition wave 
which propagates outward. This wave terminates at the point where the disk is stable. In the 
next stage of the evolution the mass in the stable region of the disk diffuses inwards, and the 
surface density of the inner region begin to increase, as does the temperature. Eventually, 
the instability is triggered anew. 

The unstable behavior is found to lie in a range of mass accretion rates, M cl ^ M ^ M c2 , 
where M cl and M c2 depends on a and a max for a given n. For larger a and/or smaller 
a max , the difference between M c2 and M cl decrease. For example, for the model parameters 
adopted (see Table 1), M cl and M c2 would be 0.11 and 0.47 M Edd for a neutron star and 
0.045 and 0.38 M Edd for a black hole of 5 M . 

Within the framework of the disk instability model, it is a general property of these 
instabilities that the frequency of the oscillation is a function of the mass accretion rate, 
with the frequency increasing for lower source intensity (assuming a direct relation between 
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intensity and mass accretion rate). Such a behavior appears to have been observed in the 
Rapid Burster (MXB 1730-335) by Lubin et al. (1992) in the hump immediately following 
the post dip phase of some Type II bursts. In addition, Lubin et al. (1992) find that the 
amplitude of the oscillations ( ^ 60%) decreased as the persistent flux declined. This is also 
consistent with the general trends exhibited by the numerical results in the previous section. 
Furthermore, the observations suggest that the level of persistent emission in the Rapid 
Burster corresponds very closely to the minimum accretion rate necessary for instability, 
M c i, since the oscillations disappear as the persistent emission decreases. If we identify the 
persistent luminosity level to correspond to M ~ 0.14M Edd , the low frequency (0.04 Hz) and 
large amplitude of luminosity oscillations imply a large n (~ 1.6) and constrain a to be 
~ 30. This value for n is similar to that inferred for a viscosity based on magnetic dynamo 
involving non-uniform rotation and cyclonic convection (Meyer & Meyer- Hofmeister 1983). 
The limit on a (i.e. a max ) may be constrained by the absence of these oscillations during 
the bump phase observed on the decline from a Type II burst (see Lubin et al. 1993). The 
intensity ratio of the bump phase to the persistent phase limits the range of mass accretion 
rates to ^ 2.5. This suggests that a max ^ 0.2. 

As shown in §4 the properties of the calculated disk oscillations are sensitive to the form 
of the viscosity parameters. Hence, the observations of oscillations from other galactic X-ray 
binary sources can provide additional constraints on the unknown parameters. Low frequency 
oscillations ~ 0.04 Hz have also been observed from the black hole candidate sources Cyg 
X-l (Angelini & White 1992; Vikhlinin et al. 1992a, 1994; Kouveliotou et al. 1992a) and 
GRO J0422+32 (Vikhlinin et al. 1992b; Kouveliotou et al. 1992b; Pietsch et al. 1993), which 
suggests that n ~ 1.2 — 1.3 and a ~ 30 for a 5M black hole. These lower values for n 
are close to those expected for a viscosity based on internal wave driven dynamo (Vishniac 
& Diamond 1992). Further long term observations, such as those carried out for the Rapid 
Burster, will be required to provide additional constraints on these parameters as well as on 

Finally, we remark that within the framework of the model proposed here one would 
expect that other X-ray binaries would also exhibit low frequency QPO phenomenon provided 
that the presence of a magnetosphere about a neutron star does not exclude the presence of 
the unstable region in the Keplerian disk. In those systems where the accretion disk extends 
to the compact object the instability may be restricted to a small range of mass accretion 
rates (as in the Rapid Burster), thus making it difficult to determine whether a given source 
can exhibit such phenomenon. In this case, transient sources would be ideal candidates for 
the search of such QPOs since the level of intensity varies over a wide range. 

The authors are grateful to A. Vikhlinin for sending a copy of his paper before publication. 
This research has been supported in part by NASA under grant NAGW-2526. 
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Table 1 

Model Sequences 



Sequence 


M(M ) 


M 


n 


a 


^max 


/ (Hz) 


AL 
17 


1 


1.4 


0.30 


1.1 


20 


1.0 




0.0 


2 


1.4 


0.27 


1.3 


20 


1.0 


0.143 


0.22 


3 


1.4 


0.30 


1.3 


20 


1.0 


0.125 


0.31 


4 


1.4 


0.21 


1.3 


20 


1.0 




0.0 


5 


1.4 


0.14 


1.6 


30 


0.2 


0.046 


0.67 


6 


1.4 


0.12 


1.6 


30 


0.2 


0.093 


0.20 


7 


1.4 


0.125 


1.6 


30 


0.2 


0.037 


0.27 














0.037 


0.12 


8 


1.4 


0.11 


1.6 


30 


0.2 




0.0 


9 


1.4 


0.47 


1.6 


30 


0.2 




0.0 


10 


10.0 


0.076 


1.2 


30 


1.0 


0.026 


0.083 


11 


10.0 


0.10 


1.2 


30 


1.0 


0.019 


0.31 


12 


10.0 


0.07 


1.2 


30 


1.0 




0.0 


13 


10.0 


0.20 


1.2 


30 


1.0 




0.0 


14 


10.0 


0.076 


1.2 


50 


1.0 


0.035 


0.056 


15 


10.0 


0.085 


1.2 


50 


1.0 




0.0 


16 


10.0 


0.071 


1.2 


50 


1.0 




0.0 


17 


10.0 


0.005-2.0 


1.1 


20 


1.0 




0.0 


18 


5.0 


0.10 


1.2 


30 


1.0 


0.05 


0.12 


19 


5.0 


0.06 


1.3 


30 


1.0 


0.036 


0.35 


20 


5.0 


0.045 


1.3 


30 


1.0 




0.0 


21 


5.0 


0.38 


1.3 


30 


1.0 




0.0 
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FIGURE CAPTIONS 



Figure 1. Curves of fi (a) and f 2 (b) used to tabulate the cooling and heating rates. The 
dotted, solid and dashed lines, correspond to the thermal equilibrium departure parameter, 
5, of 0.8, 1.0, 1.5 respectively. In this example, the central object is a 10M Q black hole, the 
viscosity parameters (n, a , a max ) are (1.2, 30, 1) and the radius is fixed at 10.57 r g . 

Figure 2. The path of the temperature and the surface density of model sequence 1 (with 
initial steady state mass accretion rate M = 0.3M Edd ) at four specified radii (in unit of 
r g ): 5.53(a), 11.77(b), 31.10(c) and 59.43(d). The empty circular points represent the initial 
steady state model and the cross points correspond to the last evolution model. The S-shaped 
curves (dashed lines) are the steady state solutions. 

Figure 3. The time variations of the bolometric disk luminosity in terms of the steady state 
value of model sequence 1 with initial steady state mass accretion rate M = 0.3MEdd- The 
disk reachs a stable state after a initial transient stage of less than 10 seconds. 

Figure 4. Same as Figure 2 but for model sequence 2. The evolution paths of the tem- 
perature and the surface density are small loops at the inner two radii and disappear at the 
other two larger radii. 

Figure 5. The time variations of the disk luminosity in terms of the steady state value. 
(a)-(c) represent model sequences 2-4 respectively. Note that sequence 2 has a weak periodic 
oscillations and sequence 4 is stable due to a lower mass accretion rate. 

Figure 6. Same as Figure 4 but for model sequence 5. In this model, the evolution paths of 
the temperature and the surface density at the inner two radii are loops, but they are much 
larger in comparision with that of sequence 2, indicating a stronger instability. Note that, 
at radius 31.10 r g (c), a small loop is present even though the initial model is located at the 
lower locally stable branch. 

Figure 7. The time variations of the disk luminosity in terms of the steady state value. 
(a)-(c) represent model sequences 5, 6 and 7 respectively. 

Figure 8. The time variations of the disk luminosity in terms of the steady state value. 
(a)-(c) represent model sequences 10, 14 and 19 respectively. 

Figure 9. The time variations of the radial distribution of the mass accretion rate of model 
sequence 19. Note that the variation is global with the same time scale at different radii, 
but the unstable region is restricted to less than 50 r g . 



18 



